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Abstract 


Development of dedicated bioenergy crop production systems will require accurate yield estimates, which will be 
important for determining many of the associated environmental and economic impacts of their production. 
Shrub willow (Salix spp) is being promoted in areas of the USA and Canada due to its adaption to cool climates 
and wide genetic diversity available for breeding improvement. Willow breeding in North America is in an early 
stage, and selection of elite genotypes for commercialization will require testing across broad geographic regions 
to gain an understanding of how shrub willow interacts with the environment. We analyzed a dataset of first-rota- 
tion shrub willow yields of 16 genotypes across 10 trial environments in the USA and Canada for genotype-by- 
environment interactions using the additive main effects and multiplicative interactions (AMMI) model. Mean 
genotype yields ranged from 5.22 to 8.58 oven-dry Mg ha“! yr _'. Analysis of the main effect of genotype showed 
that one round of breeding improved yields by as much as 20% over check cultivars and that triploid hybrids, 
most notably Salix viminalis x S. miyabeana, exhibited superior yields. We also found important variability in 
genotypic response to environments, which suggests specific adaptability could be exploited among 16 genotypes 
for yield gains. Strong positive correlations were found between environment main effects and AMMI parameters 
and growing environment temperatures. These findings demonstrate yield improvements are possible in one gen- 
eration and will be important for developing cultivar recommendations and for future breeding efforts. 
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must be established. Development and delivery of high- 
yielding, well-adapted crops is a key underlying 
assumption in the estimation of bioenergy production 
capacity in the USA (U. S. Department of Energy, 2011). 
Crop management and breeding will play crucial roles 
in meeting the challenges of producing more biomass 


Introduction 


If dedicated bioenergy crops are to play a significant 
role in climate change mitigation strategies, a clear 
understanding of yield potential on a regional basis 
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on limited land in a sustainable manner (Karp & Shield, 
2008). Shrub willow (Salix spp.) grown in short rotation 
has shown promise as a viable, regionally based feed- 
stock for marginal land due to its adaptability to cool, 
moist climates with short growing seasons (Kuzovkina 
& Quigley, 2005) and large potential for genetic 
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improvement through breeding (Smart etal., 2005; 
Smart & Cameron, 2008). This substantial yet mostly 
underexploited genetic variability among taxa is 
expected to provide a basis for developing key traits 
desirable for sustainable bioenergy production, both 
through traditional breeding and with advanced molec- 
ular techniques (Karp et al., 2011). 

Willow breeding for biomass production has been 
most extensively researched in the UK (Lindegaard & 
Barker, 1997) and Sweden (Larsson, 1998). A thorough 
account of the global breeding history of shrub willow 
is provided by Kuzovkina et al. (2008). In North Amer- 
ica, shrub willow breeding efforts began at the Univer- 
sity of Toronto in the early 1980s with collection of 
native species and exchange of plant material with the 
UK and European countries (Zsuffa, 1990), but efforts 
were focused largely on hybridizations between North 
American native species (Mosseler, 1990). In the USA, 
acquisition of plant material from collaborators in 
Canada, as well as China, Japan, New Zealand, Ukraine 
and Sweden provided the basis of a breeding program 
focused largely on novel interspecific hybridizations 
displaying heterosis for biomass yield traits (Kopp et al., 
2001; Smart et al., 2005; Smart & Cameron, 2012). Poly- 
ploidy is common in Salix and novel triploid species 
hybrids have been linked to improved yields and bio- 
mass quality in early selection trials (Smart et al., 2008; 
Serapiglia et al, 2014, 2015). Cultivar development, 
however, is a multistage screening process requiring 
substantial investments in time and resources prior to 
commercialization (Hanley & Karp, 2014). An accurate 
evaluation of cultivar yield potentials is ultimately 
assessed through testing on multiple sites with a diverse 
range of environmental characteristics. 

Biomass yield is an important factor determining the 
environmental and economic impacts associated with 
growing shrub willow. Life cycle analyses have shown 
that yield assumptions can significantly impact net 
energy ratios and greenhouse gas balances (Heller et al., 
2003; Keoleian & Volk, 2005; Caputo et al., 2014), as well 
as economic returns on investment and production costs 
(Buchholz & Volk, 2011; Hauk et al., 2014). There have 
been numerous research trials conducted over the past 
two decades in regions of the northeastern USA and 
Central and Eastern Canada that have quantified yields. 
Some studies have analyzed the differential response of 
genotype to contrasting environmental conditions 
(Labrecque & Teodorescu, 2003; Wang & Macfarlane, 
2012; Serapiglia et al., 2013; Mosseler et al., 2014a); how- 
ever, the limited number of test sites and use of diverse 
cultivars of various levels of genetic improvement 
makes it difficult to generalize specific genotypic 
responses to larger growing regions. Others have sum- 
marized mean yields from multiple test sites across geo- 


graphical regions (Kiernan et al., 2003; Volk et al., 2011), 
but due to unequal representation of genotypes among 
trials, little insight into the genotypic contribution to 
yield variability can be gained. Furthermore, efforts to 
model yields across broad geographic ranges may use 
general yield estimates from obsolete cultivars or ones 
that may not be well adapted for a particular region 
(Walsh et al., 2003; Wang et al., 2015). Plant breeders 
and agronomists are therefore challenged with assessing 
genotypic sensitivity to certain edaphic and climatic 
conditions. This process is also important for assuring 
continued improvement within a breeding program and 
for providing a basis for recommending cultivars for 
broadscale production. 

These recommendations become complicated when 
significant crossovers occur in genotype (GEN) yield 
rankings, as a response to contrasting environments 
(ENV). These genotype-by-environment (GxE) interac- 
tions are prominent and important phenomena in agri- 
culture, which present both challenges and 
opportunities for plant breeders and agronomists. Selec- 
tion and deployment of elite cultivars must be based on 
results of rigorous testing through coordinated multilo- 
cation trials, combined with appropriate statistical analy- 
ses for assessing the adaptability of genotypes and 
predicting performance in untested ENV (Annicchiarico, 
2002). A long-standing theme in plant breeding is to 
focus the search for GEN that exhibit stable yields or 
broad adaptability across ENV in the targeted growing 
region. This concept was popularized by Finlay & 
Wilkinson (1963) and Eberhart & Russell (1966), who 
developed regression parameters that seek to identify 
superior yielding GEN that maintain stable performance 
across ENV. Selection on the basis of stability can help to 
minimize the complicating effects of GxE interactions, 
adding efficiency to the selection process by focusing 
resources on material that has the best promise for wide- 
spread optimal performance (Eberhart & Russell, 1966). 
Selection based on stability should also guard against a 
potentially detrimental tendency to select GEN based on 
greater yields in only favorable ENV (Simmonds, 1991; 
Annicchiarico, 2002). However, when the GxE compo- 
nent is strong and meaningful, there is also a counter 
argument that contrasting performance among GEN can 
be capitalized upon, and breeding efforts should focus 
on specific adaptation (Cooper & Hammer, 1996; Piepho 
& Möhring, 2005). Thus, GxE interactions are viewed as 
a useful and informative aspect of cultivar testing, and 
subdividing a growing region into smaller areas and tar- 
geting GEN at those areas can improve overall yields 
(Gauch & Zobel, 1997). This argument is reinforced by 
the fact that much of the world’s crop production occurs 
on land that is less favorable than that where the crops 
were developed (Simmonds, 1991; Gauch & Zobel, 
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1997). As dedicated bioenergy crops are presumably tar- 
geted for marginal lands (Richards et al., 2014; Stoof 
et al., 2015), perhaps this concept is particularly relevant. 
The relative merits of these two perspectives, exploiting 
only broad adaptations or else both broad and specific 
adaptations, vary from case to case and depend substan- 
tially on the relative magnitudes of genotype main 
effects and GxE interaction effects. 

We present an analysis of GxE interactions in first- 
rotation yields across a network of shrub willow yield 
trials in North America covering 16 GEN and 10 ENV. 
The objectives of this study were (1) to identify shrub 
willow genotypes with broad adaptability in biomass 
yield across target growing regions in North America, 
(2) analyze GxE interactions for identifying and charac- 
terizing specific adaptation of certain genotypes and (3) 
identify edaphic and climatic variables that are most 
closely associated with GxE interaction patterns. This 
represents the most comprehensive analysis of North 
American shrub willow yields to date and will serve as 
a basis for making cultivar-site matching recommenda- 
tions and will inform future breeding efforts. 


Materials and methods 


Breeding material and yield trial network 


Foundational breeding material used in developing the GEN 
tested in these yield trials was obtained from the University of 
Toronto and from accessions of native and naturalized species 


collected in the northeastern USA and eastern Canada in the 
1980s to mid-1990s (Kopp et al., 2001; Smart & Cameron, 2008). 
Crosses were performed between 1998 and 1999 at SUNY-ESF, 
and after initial family screening in field trials for biomass yield 
traits, a group of genetically diverse individuals were deployed 
in regional yield trials. Between 2005 and 2011, 23 trials were 
established mainly in the northeastern and Midwestern USA 
and contained between six and 30 genotypes. To provide an 
unbiased comparison of genotype yields, we restricted our anal- 
ysis to 16 cultivars (Table 1) that were all present in each of 10 
environments (Table 2). The yield trials were planted between 
2006 and 2009 and hosted by institutions located in six US states 
and two Canadian provinces, and the cultivars have been placed 
into diversity groups based on pedigree. Trials were established 
and maintained generally in a consistent manner across sites 
according to a standardized protocol and methods followed 
those described in Serapiglia et al. (2013) and Volk et al. (2011). 
Conventional tillage was used to prepare the sites in either the 
fall or spring prior to planting, which generally occurred 
between May and June. All trials were planted by hand using 
dormant 25 cm cuttings sourced from nursery beds at the Tully 
Genetics Field Station of SUNY-ESF in Tully, NY. Trials were 
laid out in a double-row configuration with 1.52 m between 
double rows, 0.76 m within the double rows and 0.61 m 
between plants along the row, for a planting density of 
14 400 plants ha~'. Within each yield trial, genotype was the 
experimental treatment and the experimental units were plots 
consisting of three double rows, each 13 plants long, with the 
outer double rows serving as border rows. Genotypes were 
replicated four times in a randomized complete block design. 
Preemergence herbicides, generally oxyfluorfen (1.1 kg ai ha~') 
and simazine (2.2 kg ai ha~'), were applied prior to budbreak, 
except at Boisbriand, QC, where no herbicide was used. Periodic 


Table 1 Description of the 16 cultivars included in the genotype x environment interactions analysis of first-rotation shrub willow 


yields 

Clone ID Epithet Species / pedigree Mother Father Diversity group* Sex Ploidy+ Source 

99239-015 ‘Allegany’ S. koriyanagi x S. purpurea SH3 95058 6b F 2X Bred 

9970-036 ‘Canastota’ S. miyabeana SX61 SX64 5 M 4X Bred 

99202-004 ‘Fabius’ S. viminalis x S. miyabeana SV2 SX67 8 F 3X Bred 

9882-34 ‘Fish Creek’ S. purpurea 94006 94001 6a M 2x Bred 

99217-015 ‘Millbrook’ S. purpurea x S. miyabeana 95026 SX64 9 F 3X Bred 

9980-005 ‘Oneida’ S. purpurea x S. miyabeana 94006 SX67 9 M 3X Bred 

99113-012 ‘Onondaga’ S. koriyanagi x S. purpurea SH3 94002 6b M 2X Bred 

99201-007 ʻOtisco’ S. viminalis x S. miyabeana SV2 SX64 8 F 3X Bred 

99207-018 ‘Owasco’ S. viminalis x S. miyabeana SV7 SX64 8 F 3X Bred 

S25 ‘825’ S. eriocephala 4 F 2X Bred 

9871-31 ‘Sherburne’ S. miyabeana SX61 SX67 5 F 4X Bred 

SV1 ‘SV1’ S. x dasyclados 1 F 2X Unknown} 

SX61 ‘SX61’ S. miyabeana 5 F 4X Natural accession 
SX64 ‘SX64’ S. miyabeana 5 M 4X Natural accession 
99207-020 ‘Truxton’ S. viminalis x S. miyabeana SV7 SX64 8 M 3X Bred 

99202-011 ‘Tully Champion’ S. viminalis x S. miyabeana SV2 SX67 8 F 3X Bred 


Lower case letters indicates that these are subgroups of diversity group 6. 
*Diversity group refers to Species/ Pedigree. 


+Ploidy level estimated by flow cytometry (Serapiglia et al., 2015). 


{Collected in Ontario Canada, but possibly the cultivated hybrid S. viminalis x (S. caprea x S. cinerea) (Stott, 1991). 
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Table 2 Yield trial location characteristics. Precipitation, temperature and solar radiation data are means across four years of first rotation 


Water 
table 


Mean Solar 


Annual 


Annual 
precip. 


(mm) 


Soil 


pH 


SOM 


radiation 
(%) 


Annual 


DD 


(base 10 °C) 


G 


Year 


depth (cm) 


(MJ m7! day~') 


Tempmin (CC) 


Elevation (m) 


planted Latitude Longitude 


Code 


Location 


200 


7.07 
5.66 
6.47 
6.10 
5.25 
6.09 
6.70 
5.40 
4.80 
6.12 


4.5 


4329 


—4.1 
—1.7 


510 408 767 
1457 812 


513 


—106.61 


52.13 


2007 
2006 
2009 
2008 


Sask 


Saskatoon, SK 


45 


8.2 
3.6 


2.8 


4606 


—75.53 
—87.25 
—87.20 


43.56 


Cons 


Constableville, NY 


Skandia, MI 


30 
200 


4909 
5147 
4980 


0.1 
—0.1 


870 
1016 
1021 
1162 
1419 


822 
714 
791 
1038 
1138 


287 
222 
200 


46.36 


Skan 


45.77 


46.40 


Esca 


Escanaba, MI 
Brimley, MI 


15 


4.0 


T7 
1:2 


2.0 


—84.47 
—73.89 
—73.20 
—93.54 
—79.29 
—72.23 


2009 


Brim 


4.0 


4612 


30 
114 


45.63 
44.01 


Boisbriand, QC Bois 2007 


Middlebury, VT 
Waseca, MN 


30 
30 
30 
45 


6.8 


4612 


2007 
2006 
2008 


Midd 


5.7 
3.6 
3.4 


4944 
5082 
5171 


1,7 
4.9 


1459 
1477 
1487 


843 
909 
1274 


349 
255 


44.06 


Wase 


42.44 


41.80 


Fred 
Stor 


Fredonia, NY 
Storrs, CT 


5.4 


198 


2009 


GDD, growing degree-days; SOM, soil organic matter. 


mechanical or spot chemical weed control was performed as 
needed. After the first year of growth, all stems were cut back 
during dormancy to promote coppice regrowth the following 
spring, at which time 112 kg N ha™' was applied, normally as 
ammonium sulfate, except for the trial located in Saskatoon, SK, 
where no fertilizer was applied. Harvests were conducted three 
years postcoppice (4-year-old root systems) during dormancy, 
except at Brimley, MI, which was harvested in July at 3.5 years 
postcoppice due to extreme moisture conditions in the field. 
Two to four plants on each end of the middle double row of each 
plot were excluded from harvest measurements to avoid plot-to- 
plot edge effects, resulting in 18-22 plants available for measure- 
ments. All stems from the 18-22 plants from the middle double 
row of each plot were cut with brush saws or cut and chipped 
with a mechanical harvester and weighed in a bin with weigh 
cells in the field. A subsample from each plot consisting of either 
whole stems or chips was collected, weighed fresh, dried at 
65 °C to a constant weight and reweighed to determine moisture 
content. Moisture content was used to calculate dry matter yield 
for each plot based on the area occupied by the harvested plants 
across a 3-year rotation. All yields reported here are expressed 
as oven-dried Mg ha™! yr~'. Survival data was also recorded on 
the harvested plants at the time of harvest, except Boisbriand, 
QC, where no survival data was collected, but survival was 
>90% for each plot (M. Labrecque, Personal observation) and 
Savoy, IL where survival was assessed in the middle of the sec- 
ond rotation, but most of the mortality occurred early in the first 
rotation (G. Kling, personal observation). 


Site environmental characteristics 


Daily temperature and precipitation data for all 4 years of the 
harvest cycle were obtained from weather stations nearest to 
each trial location with the most complete records using publi- 
cally accessible databases, National Oceanic and Atmospheric 
Agency, National Centers for Environmental Information 
(NOAA NCEI, 2015) for US trials and Canadian National Cli- 
mate Data (Environment Canada, 2015) for CA trials. Daily 
solar radiation estimates at a 1° by 1° grid scale were obtained 
from National Aeronautics and Space Adminstration, Predic- 
tion of Worldwide Energy Resource (NASA POWER, 2015). 
Soil samples were generally collected at the time of planting, or 
occasionally soon after harvest. Due to some differences 
between soil extraction methods, only pH (1 : 1 soil/water by 
weight) and % organic matter (typically by loss on ignition) 
were considered for data analysis. 


Dataset refinement 


Initially, our dataset was comprised of 640 independent obser- 
vations where the harvested plot was the experimental unit. 
While overall mean survival was >90% across the yield trial 
network, some trials and individual plots experienced greater 
mortality. Two trials, Escanaba, MI, and Saskatoon, SK, had 
mean survival values below 80%. Damage from herbicide drift 
appeared to be the main cause of mortality in Saskatchewan 
(Amichev et al., 2015), while at Escanaba mortality could not be 
associated with any particular issue (R. Miller, personal 


© 2016 The Authors. Global Change Biology Bioenergy Published by John Wiley & Sons Ltd., doi: 10.1111/gcbb.12344 


GxE INTERACTIONS IN SHRUB WILLOW HYBRIDS 5 


observation). In total, 39 independent observations (experimen- 
tal units) with >65% mortality at the time of harvest were 
removed from the original dataset of 640 observations. Data 
were also inspected for extreme moisture content values, which 
can impact dry matter yield calculations. Two likely sources of 
error were 1) premature removal of samples from the drying 
oven, resulting in excessively low moisture content estimates, 
and 2) extraneous moisture contamination at harvest during wet 
conditions, leading to excessively high moisture values. Mois- 
ture values that were greater or less than three standard devia- 
tions of the mean across all samples were considered outliers 
and were removed from the dataset. To retain a yield value for 
the particular observation, the mean of the remaining three 
replicates was applied to the plot fresh weight of the outlying 
moisture content. This was performed for 15 samples in total 
from four of nine trials, which represented 2.3% of the total 
observations considered in this analysis. For the Boisbriand, QC 
trial, a single wood sample was collected for each cultivar and 
the moisture content was applied to all four replicates of that 
cultivar. The mean moisture content of these samples was 
reported at 32.7% (SD 2.90%), and it was assumed that all sam- 
ples had not dried sufficiently, considering the overall mean 
moisture content across all other trials was 46.9% (SD 4.77%). A 
correction factor equaling the difference between these two 
mean values was added to each of the observations originally 
reported for the Boisbriand trial. This value was 14.2%, which 
was also very close to the value of 3 standard deviations of the 
overall mean (14.3%). After accounting for survival and mois- 
ture content adjustments, 601 observations were available for 
analysis from the original 640 from 16 GEN in 10 ENV. 


Statistical analysis 


We chose to analyze GxE interactions in our first-rotation yield 
dataset using the additive main effects and multiplicative inter- 
actions (AMMI) model. The AMMI model is a combination of 
analysis of variance (anova) and principal components analysis 
(PCA), where the GxE interaction, contained in the residual of 
the additive GEN and ENV main effects model, is subjected to 
PCA and the interaction sum of squares (SS) are partitioned into 
a series of interaction principal component (IPC) axes, where 
IPC1 is the first interaction PC axis, and so on. The AMMIO 
model indicates that no IPC axes are included and the model is 
equivalent to the additive (main effects only) anova. AMMI1 
includes the first IPC axis, and so on, while AMMIF is the full 
model that includes all axes and is equal to the raw data. The 
maximum number of IPC axes for a given dataset is one less the 
minimum number of GEN or ENV (Gauch, 1992; p. 85), but typi- 
cally the majority of the interaction signal is captured in the first 
few IPC axes, with later IPC axes containing increasing amounts 
of interaction noise, and decreasing amounts of signal (Gauch, 
2012). Therefore, a more parsimonious model containing a lower 
number of axes is favored and the remaining axes are relegated 
to a pooled residual. The statistical significance of each axis is 
often assessed with an F-test where the degrees of freedom (DF) 
for each axis are calculated according to Gollob (1968) (see also 
Gauch, 1988). More recently, Piepho (1995) demonstrated that 
the traditional F-test can be too liberal and suggested a more 


conservative Fr-test for determining the significance of each IPC 
axis. Also, Gauch (1992, p. 147) suggested a simple test for 
model diagnosis and selection, which involves estimation of the 
GxE noise SS by multiplying the GxE DF by the mean square 
error. Consequently, the signal GxE SS can be estimated by sub- 
tracting the GxE noise SS from the GxE total SS. Therefore, IPC 
axes can be assessed by the amount of GxE signal SS that are 
recovered, instead of the total GxE SS. Selection of higher order 
models that include greater numbers of axes must be weighed 
not only in terms of statistical significance, but also in terms of 
practicality, because higher order AMMI models and especially 
the noisy AMMIF tend to have a large roster of genotypes that 
win in at least one environment, and hence, such models pro- 
duce an unmanageable number of mega-environments (Gauch, 
2013). Parsimonious models are often preferred, and AMMI 
models are most useful when the multiplicative terms have agri- 
cultural interpretability (Gauch, 2013). Equation 1 gives the gen- 
eral form of the AMMI model: 


Y ger = pr Ag + Be + Ln An? gnen T Pree Ar Eger (1) 


where Yger is the yield of the gth genotype in the eth environ- 
ment for the rth replicate, u is the grand mean, a, is the geno- 
type g mean deviation (genotype mean minus grand mean), pe 
is the environment e mean deviation, 2; is the singular value for 
nth IPCA axis, ygn is the genotype g eigenvector value for IPCA 
axis n, den is the environment e eigenvector value for IPCA axis 
n, Pge is the residual, and éger is the experimental error. 

The AMMI analysis was performed using MATMODEL 
V3.0, an open source statistical program designed specifically 
for the analysis of GxE interactions (Gauch & Furnas, 1991; 
Gauch, 2007). MATMODEL performs the combined ANova/ 
PCA and also delineates mega-environments, which allows for 
the exploration of specific adaptation to particular edaphic or 
climatic conditions (Gauch & Zobel, 1997). 

MATMODEL does not analyze the experimental design; 
hence, we first analyzed the experimental design using PROC 
MIXED in SAS® version 9.4 (SAS Institute Inc., 2013), with 
block nested within ENV, and coded as a random effect. The 
main effects of GEN and ENV and the GxE interaction term 
were considered fixed effects. The TYPE3 option in PROC 
MIXED was invoked to obtain expected SS for the fixed effects. 
Tukey’s studentized range (HSD) post hoc test was performed 
for means separation among GEN. The least squared means 
from the SAS output were then supplied to MATMODEL to 
perform the PCA of the interaction and for the mega-environ- 
ment analysis. The random error variance from the SAS PROC 
MIXED output was used to calculate the F-tests in the AMMI 
analysis. Finally, the environment main effects and IPC scores 
were used in correlation analyses with the various environmen- 
tal variables using PROC CORR in SAS® version 9.4. 


Results 


ANOVA and genotype main effects 


The grand mean for first-rotation yields for this base 
case dataset based on the mixed model analysis of vari- 
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ance was 6.96 Mg ha_' yr '. The main effects of ENV 
and GEN and the GxE interaction were all highly 
significant (P < 0.0001; Table 3) and accounted for 82, 6 
and 12% of treatment SS, respectively. The blocking 
effect (nested within ENV) accounted for a relatively 
large amount of random error variance (Table 3). 
Although the main effect of genotype accounted for 
only 6% of the treatment SS, it contains important infor- 
mation about patterns of broad adaptation. The cultivar 
Salix viminalis x S. miyabeana ‘Fabius’, a triploid hybrid 
in diversity group 8 (DG8) was the overall greatest 
yielding GEN across the 10 ENV tested, with a mean 
yield of 8.58 Mg ha ' yr ' (Table 4). Two other Salix 
viminalis x S. miyabeana triploid hybrids in DG8, 
‘Otisco’ and ‘Tully Champion’, ranked 2nd and 4th in 
overall yield, respectively (Table 4, Fig. 1a). Two of the 
top five GEN were tetraploid S. miyabeana (DG5), 
including ‘Canastota’ a progeny selected from a cross 
between ‘SX61’ and ‘SX64’. Two triploid S. pur- 
purea x S. miyabeana cultivars from DG9, ‘Oneida’ and 
‘Millbrook’, were ranked 6th and 7th overall, respec- 
tively (Fig. 1a). Triploids showed the greatest relative 
yields (GEN yield/ENV mean yield averaged over all 
ENV), with the exception of two cultivars, ‘Truxton’ 
and ‘Owasco’, which performed at or below ENV 
means (Fig. 1b). The cultivars S. miyabeana ‘SX61’, 
‘SX64’ and S. x dasyclados ‘SV1’ were cultivars that 
showed promise in earlier trials, and subsequently were 
used as check cultivars in this yield trial network. The 
top two improved cultivars, ‘Fabius’ and ‘Otisco’, per- 


Table 4 Adjusted mean first-rotation shrub willow yields (Mg ha 


= 


Table 3 Mixed model analysis of variance for first-rotation 
yields of 16 shrub willow genotypes in 10 environments show- 
ing F-test results for fixed effects and variance components for 
random effects 


Source DF SS MS F Value Pr>F 
Fixed 
ENV 9 5247.57 583.06 28.01 <0.0001 
GEN 15 381.03 25.40 9.19 <0.0001 
GxE 135 765.39 5.67 2.05 <0.0001 
DF VC Pct 
Random 
BLK (ENV) 30 1.27 31.43 
Error 411 2.76 68.57 


DF, degrees of freedom; SS, sum of squares; MS, mean square; 
ENV, environment; GEN, genotype; BLK, block; VC, variance 
component; Pct, percent of total variance. 


formed better than all three check cultivars, and four 
other improved cultivars performed better than the 
check clone mean yield of 7.14 Mg ha! yr. The GEN 
with the lowest mean yields were diploid S. eriocephala 
‘$25’ (DG4), a North American native willow species, 
and two hybrid cultivars of S. koriyanagi x S. purpurea 
‘Onondaga’ and ‘Allegany’ (DG6b). ENV mean yields 
ranged from 2.57 to 11.3 Mg ha! yr, with the great- 
est yields in eastern USA and Canada, and the lowest 
yields occurring in the Upper Peninsula of MI, USA, 
and Saskatoon, SK, CA. 


yr_') for 16 genotypes in 10 environments in North America 


Environments* 

Epithet Midd Bois Stor Wase Fred Cons Esca Sask Skan Brim MEAN 
‘Fabius’ 13.96 13.17 13.27 10.01 8.19 8.76 8.83 4.38 1.95 3.25 8.58 
‘Otisco’ 11.95 12.39 11.10 7.66 7.54 6.55 9.38 4.46 3.18 2.91 7.71 
‘SX64’t 11.11 13.61 9.15 10.13 8.01 8.61 5.94 3.34 3.21 3.08 7.62 
‘Tully Champion’ 13.34 13.15 6.65 7.57 8.22 5.61 7.51 5.80 3.75 3.70 7.53 
‘Canastota’ 14.56 13.36 8.78 8.98 7.71 7.75 5.20 3.53 2.29 2.47 7.46 
‘Oneida’ 10.63 8.85 12.31 9.80 6.50 8.89 6.91 4.67 3.25 2.41 7.42 
‘Millbrook’ 10.77 10.04 9.23 9.11 8.00 8.11 7.96 4.13 2.24 2.63 7.22 
‘SX61’ 11.97 10.45 9.44 9.62 6.78 7.05 7.89 2.75 3.05 3.14 7.21 
‘Fish Creek’ 11.18 10.78 9.77 9.07 7.59 8.17 6.30 4.03 2.00 1.74 7.06 
‘Truxton’ 12.32 9.95 9.84 8.96 8.32 6.59 5.71 3.60 3.40 1.74 7.04 
‘Sherburne’ 11.92 10.72 9.82 8.61 6.53 7.02 5.50 3.32 1.89 2.30 6.76 
‘SVY 10.92 11.54 6.03 4.11 7.62 7.81 9.02 3.28 2.47 3.23 6.60 
‘Owasco’ 10.13 8.75 7.18 8.85 8.29 6.32 6.44 4.20 3.15 2.56 6.59 
‘Allegany’ 8.31 9.31 6.88 7.51 5.26 6.32 6.94 3:10 2.21 2.13 5.80 
‘Onondaga’ 8.36 7.54 5.54 7.34 5.10 6.84 6.67 3.60 2.06 2.02 5.51 
‘$25’ 9.44 8.52 5.39 6.60 6.52 3.87 3.90 3.81 2.39 1.74 5.22 
MEAN 11.30 10.76 8.77 8.37 7.26 7.14 6.88 3.87 2.66 2.57 


*Environment names are truncated to the first four letters, with full environment names given in Table 2. 


+Check cultivars are underlined. 
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Fig. 1 First-rotation shrub willow mean yields (a) and yields 
as a proportion of environment means (b) for 16 genotypes in 
10 environments in North America. Bars are shaded according 
to the diversity groups (pedigrees) presented in Table 1. Bars 
in (a) with different letters indicate significant (P < 0.05) differ- 
ences according to Tukey’s studentized range test (HSD). Error 
bars in (a) are + one standard error of the mean. Proportions in 
(b) were calculated as the yield of a genotype divided by the 
environment mean, averaged across the 10 environments. 


AMMI analysis 


The AMMI2 model anova revealed that the main effects 
of GEN and ENV plus the first two IPC axes accounted 
for 95.5% of the treatment SS with only 68 (of 159) treat- 
ment DF (Table 5). Using the error variance obtained 


from the SAS PROC MIXED output (Table 3), both the 
F-tests and the more conservative Fr-tests proved IPC1 
and IPC2 to be significant (P < 0.0002); however, all 
subsequent axes were not significant and were therefore 
pooled into a discarded residual. The GxE interaction 
term is inherently a mixture of true signal in the data 
G.e., systematic rank changes amount GEN) and noise. 
The amount of GxE noise can be estimated by multiply- 
ing the GxE DF by the mean square error (Gauch, 
2013). Using the mean square error variance obtained 
from the SAS PROC MIXED output (Table 3), the GxE 
noise SS for this dataset was estimated to be 373 (46%), 
while the GxE signal accounts for 438.5 (54%) of the 
total GxE SS. IPC1 captured 291 SS, or 66.4% of 
the GxE signal SS. The cumulative SS accounted for by 
the first two IPC axes is 507, which is greater than the 
estimated GxE signal SS, suggesting that IPC2 contains 
a combination of mostly signal and some noise. There- 
fore, AMMII does a reasonable job of capturing the 
majority of the GxE signal in a parsimonious model. 
However, in terms of accurately describing our data, 
AMMI2 might be slightly better, although the more 
important consideration is that AMMII (and AMMI2) is 
considerably more accurate than the raw data AMMIF 
as its additional IPC3 and higher components (which 
are combined in the residual in Table 5) capture a SS of 
304, which is mostly noise. 

A useful tool for interpreting the AMMI analysis is 
the AMMI]1 biplot, where the additive main effects of 
GEN and ENV are plotted on the x-axis in units of yield 
(Mg ha"! yr~'), and the IPC1 scores are plotted on the 
y-axis, which have units expressed as the square root of 
the yield (Gauch, 1992; p. 85). By representing both the 
main effects and the majority of the interaction signal in 
one projection, the AMMI1 biplot captures 92.3% of the 
treatment SS, and the relationships between main effects 
and interaction can be observed simultaneously (Fig. 2). 
The vertical reference line represents the grand mean 
yield of the dataset, and the horizontal line is placed at 
zero for the IPC scores, where points farther away from 
this line indicate GEN or ENV with larger interactions. 
Genotypes that occur close to the horizontal line can be 
regarded as having relatively stable yields across ENV. 
Genotypes and ENV having the same sign for their 
IPC1 scores have positive interactions, whereas opposite 
signs give negative interactions. For instance, ‘Fabius’ 
was the highest yielding GEN, and consequently, it is 
farthest along the right side of the x-axis. It also had one 
of the largest IPC1 scores, generating a positive interac- 
tion in Storrs, CT, but a negative interaction in Skandia, 
MI. Indeed, the lowest yield for ‘Fabius’ of 
1.95 Mg ha! yr! occurred at that location, in compar- 
ison with Storrs, CT, where ‘Fabius’ yielded 
13.27 Mg ha ! yr! (Table 4). The IPC1 scores for ‘SV1’ 
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Table 5 Additive main effects and multiplicative interactions (AMMI) analysis of variance for first-rotation yields of 16 shrub willow 


genotypes in 10 environments showing the first two IPC axes 


Source DF SS MS 

TRT 159 6774.19 42.60 
GEN 15 454.50 30.30 
ENV 9 5508.20 612.02 
GxE 135 811.49 6.01 
IPC1 23 290.96 12.65 
IPC2 21 216.04 10.29 
Residual 91 304.49 3.35 


F value* Pe SF Fr Value Pr> Fr 
15.42 <0.0001 
10.97 <0.0001 
221.51 <0.0001 
2.18 <0.0001 
4.58 <0.0001 2.18 <0.0001 
3.72 <0.0001 1.68 <0.0002 
1.21 0.99 1.21 0.11 


DF, degrees of freedom; SS, sum of squares; MS, mean square; TRT, treatment; ENV, environment; GEN, genotype; BLK, block. 
*The error mean square from Table 3 was used in calculations of F and Fr values. 
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Fig. 2 AMMI1 biplot for the shrub willow yield trial network. Genotype (e) and environment (A) means are on the x-axis, and the 
IPC1 scores are shown on the y-axis. Genotype and environment names are truncated to the first four letters, and full names can be 
found in Tables 1 and 2, respectively. The vertical reference line represents the grand mean, while the horizontal line crosses at an 


IPC score of zero. 


and “Tully Champion’ are similar and were the most 
extreme negative scores out of the 16 GEN, so their 
interaction patterns are the opposite of ‘Fabius’. How- 
ever, ‘Tully Champion’ has a substantially greater mean 
yield compared to ‘SV1’, and it was the top-yielding 
cultivar in many low-yielding ENV (Table 4). 

When IPC2 is significant (Table 5), the AMMI2 biplot 
can be used to investigate the interaction pattern in 
IPC1 and IPC2 together (Fig. 3). The reference lines in 
Fig. 3 are drawn through zero for both axes. The cross- 
ing of these two lines in the middle of the graph indi- 
cates no interaction, and therefore, a GEN close to this 
point would be characterized as having stable yields 
across ENV, and for this dataset, the triploid ‘Otisco’ 
and the tetraploids ‘SX61’ and ‘SX64’ are closest to the 


origin. ‘Fabius’ and ‘Tully Champion are on the oppo- 
site ends of IPC1 scores, but have very similar IPC2 
scores. In contrast, ‘Onondaga’ and ‘Canastota’ were on 
the opposite extremes of IPC2 scores, with ‘Canastota’ 
having a strong positive interaction with Boisbriand, 
QC in terms of IPC2, as their IPC scores are of the same 
sign. Unlike Fig. 2, there is no information about main 
effects in Fig. 3, only interactions, but it is useful for 
understanding the relationships between AMMI1 and 
AMMI2. 


Mega-environment delineations 


Another useful graphical representation of AMMII 
analysis is a plot of the AMMI1 nominal yields (AMMI1 
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Fig. 3 AMMI2 biplot for the shrub willow yield trial network 


. Genotype (e) and environment (A) IPC1 scores are shown on the 


x-axis, and IPC2 scores are on the y-axis. Genotype and environment names are truncated to the first four letters, and full names can 
be found in Tables 1 and 2, respectively. The vertical and horizontal reference lines are drawn through zero for both axes. 
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Fig. 4 Shrub willow genotypic responses according to the AMMI1 model. The 10 environment IPC1 scores are represented on the x- 
axis, while the genotype nominal yields (AMMI1 model yields without the environment deviations) are shown on the y-axis. All 16 
genotypes are represented individually by a straight line. Only relevant genotypes have been labeled. Check cultivars are solid black 
lines. Environment names are truncated to the first four letters, and full names can be found Table 2. 


model yields for each GEN, without the ENV deviation, 
see equation 1), against the ENV IPC1 scores (Gauch & 
Zobel, 1997). Figure 4 shows nominal yields for the 16 
GEN across the 10 ENV plotted as straight lines, and 
the vertical positions of the lines relative to one another 
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at a given ENV IPCI score indicate the yield rankings. 
This plot is very useful for visualizing GEN crossovers, 
and for ease of interpretation, only the most relevant 
GEN have been highlighted and their names provided. 
For instance, ‘Fabius’ had the greatest yields in five 
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ENV ranging in IPC1 scores from —0.38 to 2.16. But for 
the other five ENV with scores —0.50 and lower, a 
switch point occurs where ‘Tully Champion’ outper- 
forms ‘Fabius’ according to the AMMI1 model pre- 
dicted yields. This situation illustrates the concept of 
narrow adaptation, where GEN perform similarly in a 
group of ENV with similar characteristics. 

The degree of environmental sensitivity of a GEN can 
also be inferred by observing the slopes in Fig. 4. 
‘Fabius’ and ‘Tully Champion’ have steep, but opposite 
slopes due to their opposite interactions. ‘Otisco’ was 
ranked 2nd in overall yield and was relatively insensi- 
tive to ENV, as evidenced by its near zero slope (Fig. 4). 
Also, similarities among GEN can be observed where 
slopes are nearly parallel. The line just below ‘Fabius’ 
on the right side of Fig. 4 is that of ‘Oneida’, which is 
similar to ‘Fabius’, suggesting that they reacted to ENV 
similarly. However, the mean yields of ‘Fabius’ were 
considerably greater than ‘Oneida’, making ‘Fabius’ an 
obvious superior choice. Similarly, ‘SV1’ had an interac- 
tion pattern that mirrored ‘Tully Champion’, but again, 
the predicted AMMI1 yields were lower for ‘SVT. 

AMMI can also be used to delineate mega-environ- 
ments, which are defined as subdivisions of a crop’s 
potential growing range that share similar genotype 
winners, presumably due to similar biotic and abiotic 
stresses, but that are not necessarily contiguous (Gauch 
& Zobel, 1997; Yan et al., 2000). Subdivision of a crop’s 
growing region can help to improve yields by targeting 
GEN with specific adaptation to ENV where they are 
likely to perform best (Gauch & Zobel, 1997), especially 
if there are identifiable environmental or biological pat- 
terns associated with the IPC analysis that can be 


extended to ENV beyond those where crop yields have 
been tested (Gauch, 2013). MATMODEL uses key 
switch points in nominal yields among winning GEN 
across IPC1 scores to delineate mega-environments. For 
this shrub willow yield dataset, the main switch point 
for the AMMI1 model is clearly illustrated in Fig. 4, 
where ‘Fabius’ and ‘Tully Champion’ switch between 
first and second rank at a value of about —0.44 along 
the ENV IPC1 range. Therefore, AMMI1 defines two 
winners which happen to divide the 10 ENV equally 
between them, where ‘Fabius’ is declared the overall 
winner, but ‘Tully Champion’ is superior in the lower 
yielding ENV. 

As the AMMI analysis consists of a family of models, 
with higher order models incorporating greater num- 
bers of interaction components, higher order models 
generally result in an increased number of mega-envir- 
onments and winners (Gauch, 2013). To illustrate this 
point, Table 6 shows the top five rankings for three 
AMMI model family members. The ENV are listed by 
IPC1 scores. AMMI1 divides the growing region into 
two mega-environments as described above, with only 
two-first place winners, ‘Fabius’ and ‘Tully Champion’. 
‘Fabius’ takes second place in most ENV where ‘Tully 
Champion’ was ranked first. ‘Oneida’ ranks behind 
‘Fabius’ in more favorable ENV, while ‘Otisco’ ranks 
third in many of the less favorable ENV. The main dif- 
ferences between AMMI1 and AMMI? are that Bois- 
briand, QC is declared as a separate mega-environment 
with ‘Canastota’ the top-ranking GEN, and to some 
degree, ‘Fabius’ is relegated to lower ranks in the less 
favorable ENV. The AMMIF model, which is equivalent 
to the raw data, is also included for comparisons. 


Table 6 Rankings of the top five genotypes for in 10 environments (ENV) based on three AMMI models. AMMI]1 includes only the 
first IPC axis and AMMI2 includes both the first and second IPC axes which were both significant based on the Fr-test. AMMIF repre- 


sents the raw data 


AMMI1 Rank AMMI2 Rank AMMIF Rank 

ENV* 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 
Stor Fabit (Onei Fish SX64 Sher Fabi Onei Fish SX64 Sher Fabi Onei Otis Trux Sher 
Wase Fabi Onei SX64 Otis Cana Fabi Onei Mill Fish Otis SX64 Fabi Onei SX61 Mill 
Fred Fabi Onei Otis SX64 Cana Fabi Onei Otis Mill SX64 Onei Fabi SX64 Fish Mill 
Midd Fabi Otis SX64 Onei (Cana Fabi Cana Tull SX64 Otis Cana Fabi Tull Trux SX61 
Bois Fabi Tull Otis SX64 Cana (Cana Tull Fabi SX64 Otis SX64 Cana Fabi Tull Otis 
Esca Tull Fabi Otis SX64 Cana Tull Fabi Otis Cana SX64 Trux Owas Tull Fabi SX64 
Sask Tull Fabi Otis SV1 SX64 Tull Otis Fabi SV1 Mill Tull Onei Otis Fabi Owas 
Cons Tull Fabi Otis SV1 SX64 Tull Otis Mill SV1 Owas Otis SV1 Fabi Mill SX61 
Skan Tull Fabi Otis SV1 SX64 Tull Otis SV1 Mill Fabi Tull Trux Onei SX64 Otis 
Brim Tull SV1 Fabi Otis SX64 Tull SV1 Otis Fabi SX64 Tull Fabi SV1 SX61 SX64 


*ENV, Environment. Names are truncated to the first four letters. Full names can be found in Table 2. 
+Genotype names are truncated to four letters, with full names presented in Table 1. The overall top-yielding cultivar, ‘Fabius’ is 
underlined throughout the table, and the second mega-environment winner from AMMI1, ‘Tully Champion’ is italicized throughout 


the table. 
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Should the full data be accepted as most accurate, it 
would result in a very complex array of mega-environ- 
ments with seven different cultivars being declared win- 
ners in seven narrowly defined subregions. AMMI1 
with its simplified mega-environment winner pattern is 
likely best for making cultivar selections far less compli- 
cated. AMMI2 is likely the most accurate choice for 
modeling overall genotypic performance in this dataset, 
but the inclusion of one winner in one ENV complicates 
recommendations and justification for its inclusion 
would likely need to be verified with more testing. 


Edaphic and climatic variables 


The AMMI analysis of the GxE interaction operates 
solely on yield data, with an implicit interpretation that 
general patterns in yields are a reflection of underlying 
environmental characteristics of test locations. However, 
it may be desirable to associate yield patterns with 
known environmental characteristics, and IPC scores 
can be used in simple correlations when environmental 
variables for test locations are available (Van Eeuwijk, 
1995). ENV yield was significantly (P < 0.05) correlated 
with mean annual growing degree-days and marginally 
significantly (P < 0.10) associated with longitude and 
growing season precipitation (Table 7). ENV IPC1 
scores were also significantly correlated with growing 
degree-days, mean annual minimum temperatures and 
latitude (Table 7), which would covary with each other, 
but suggest a strong relationship between interaction 
patterns and temperature. IPC2 scores were positively 
correlated with elevation (Table 7). 


Expanded datasets and analyses 


As noted earlier, our yield trial network contained more 
GEN and test ENV than those included in the above 


Table 7 Pearson product-moment correlation coefficients between mean yields (Mg ha 


variables for the 10 environments in the yield trial network 


AMMI analysis (601 observations of ~1530 available). 
However, the MATMODEL program allows for imputa- 
tion of missing treatment cells in the two-way table of 
means using the expectation maximization (EM) algo- 
rithm (Gauch & Zobel, 1990). To more fully explore the 
yield database, two additional datasets were assembled 
that expanded the original complete set of 16 GEN in 10 
ENV. The first dataset had an additional 6 GEN that 
included 10 missing treatments (4.5%) of a possible 220, 
for a total of 22 GEN in 10 ENV. The second dataset 
added 6 ENV with 15 missing treatments (5.9%) of 256. 
It has recently been shown that using the EM method 
for data imputation with missing data proportions of 5— 
10% do not significantly reduce the predictive ability of 
the AMMI model, especially when a small number of 
IPC axes (one or two) are chosen (Rodrigues et al., 2011; 
Paderewski & Rodrigues, 2014). The anova tables and 
mega-environment winners for the two expanded data- 
sets are provided in Table S1. The anova for the two 
expanded datasets show greater overall variance in the 
treatments, as expected, but the GxE signal SS were 
similar at 38.7% for the increased GEN and 36.7% for 
the increased ENV scenarios (Table S1). The IPC1 axes 
accounted for 82.4 and 76% of the GxE signal SS for the 
expanded GEN and expanded ENV scenarios, respec- 
tively. Interestingly, the mega-environment analyses 
remain remarkably similar to those of the base case sce- 
nario (Table $2). Even with six additional GEN, ‘Fabius’ 
remains the winner in six ENV and ‘Tully Champion’ is 
the winner in four of the lowest yielding ENV. When 
six ENV were added for a total of 16, ‘Fabius’ wins in 
seven and ‘Tully Champion’ in eight (Table S2). ‘SVT 
wins in one environment, Belleville, NY, but ‘Tully 
Champion’ is a close second, and often combining a 
small mega-environment with a larger one occupied by 
a close winner is justified (Gauch, 2013). The addition of 
these two scenarios reinforces the superiority of ‘Fabius’ 


~! yr~!), IPC scores and edaphic and climatic 


Variable Yield IPC1 IPC2 
Latitude —0.534 —0.619* 0.098 
Longitude 0.587* 0.300 —0.503 
Elevation (m) —0.496 —0.204 0.785** 
Mean annual precipitation (mm) 0.527 0.346 —0.221 
Mean April—October precipitation (mm) 0.578* 0.525 —0.195 
Mean annul growing degree-days (base 10 °C) 0.650** 0.827** —0.287 
Mean annual minimum temperature (°C) 0.394 0.729** —0.166 
Mean annual solar radiation (MJ m7! day~') —0.038 0.479 0.238 
Soil organic matter (%) 0.261 —0.116 —0.013 
Soil pH 0.028 —0.191 —0.320 
Depth to water table (cm) —0.131 —0.300 0.009 


*P < 0.10; **P < 0.05. 
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and ‘Tully Champion’ relative to the cultivars tested in 
the analysis. 

Incidentally, MATMODEL can also provide the Fin- 
lay—Wilkinson joint regressions analysis (Zobel et al., 
1988; Gauch, 2007), which we applied to the original 
dataset of 16 GEN in 10 ENV. The joint regression anal- 
ysis only accounted for 23.2% of the GxE SS, compared 
to 35.9% for AMMI IPC1 (Table S3). Based on GEN 
slopes, the joint regression analysis of the GEN with the 
best combination of stable slopes and greater yields 
would be as ‘Otisco’ and ‘SX64’, with the high yields 
and slopes slightly above 1, and ‘Tully Champion’ and 
‘Oneida’ with somewhat lower yields and slopes 
slightly less than 1 (Table S3). These results are some- 
what in agreement with AMMI, but selection based on 
stability would likely rule out the potential gains from 
including ‘Fabius’ in optimal locations. 


Discussion 


Our analysis of this North American shrub willow yield 
trial network data has demonstrated the importance of 
GxE interactions in short-rotation woody crop produc- 
tivity. Our evaluation of first-rotation yields using the 
AMMI model has accomplished the main objectives of 
the study, which were to identify GEN with broad 
adaptation, to define subregions where particular GEN 
exhibit specific adaptation and to identify environmen- 
tal variables that help to explain patterns in yield. 
Despite the proven success of AMMI in accurately diag- 
nosing GxE interactions from a range of agronomic and 
horticultural crops, we have not seen any other pub- 
lished reports of its use in short-rotation coppice 
research. AMMI helped to resolve signal from noise, 
and we selected the more parsimonious model, AMMI1, 
as the best representation of the data, simplifying the 
mega-environment analysis with just two winning 
GEN, ‘Fabius’ and ‘Tully Champion’, as opposed to the 
raw data with seven of the 16 GEN in our dataset win- 
ning in at least one environment (Table 6). In further 
support of the AMMI1 analysis winners, inclusion of 
additional GEN or additional test locations did not 
change the mega-environment analysis outcomes for 
AMMII (Table S1). 

‘Fabius’ was the top-yielding cultivar with an overall 
mean of 8.58 Mg ha! yr“, setting a standard for opti- 
mal yields and broad adaptability. Also, because of its 
positive interaction with high-yielding ENV, it is well 
adapted to favorable growing conditions. “Tully Cham- 
pion’ had an opposite pattern of interaction and appears 
to be adapted to less desirable growing conditions. 
‘Tully Champion’ and ‘SV1’ interacted with the test 
ENV similarly (Fig. 4), but ‘Tully Champion’ consis- 
tently outperformed ‘SVT make it a much better check 


cultivar for evaluating future breeding material. Simi- 
larly, ‘Otisco’ and ‘SX64’ had fairly stable yields, show- 
ing little interaction with environment regardless of 
quality, but ‘Otisco’ was ranked second in overall mean 
yield. All three of these superior cultivars (‘Fabius’, 
‘Tully Champion’ and ‘Otisco’) are of the triploid pedi- 
gree of S. viminalis x S. miyabeana, and they outper- 
formed the mean of the check cultivars by 5 to 20%, 
demonstrating that yield gains can be achieved through 
novel hybridization in shrub willow. This affirms the 
findings of Serapiglia et al. (2014) for early selections of 
triploids on a single site, but extends it more broadly 
across a large set of diverse ENV. 

It should be noted that ‘Fabius’ and “Tully Champion’ 
are siblings from the same cross, which makes their 
divergent patterns in yield across our sites rather 
remarkable. These two were selected after an evaluation 
of the family in a replicated trial at a single location in 
central NY (Smart et al., 2008). There was considerable 
variability in growth potential across this family, and 
clearly the single test location was not adequate to 
detect differiential response to broad environmental 
conditions. The cultivar with the most stable yields in 
our dataset, ‘Otisco’, shares the same S. viminalis 
mother, ‘SV2’, with ‘Fabius’ and ‘Tully Champion’ 
(Table 1). ‘Owasco’ and ‘Truxton’ are also siblings that 
belong to the Salix viminalis x S. miyabeana pedigree 
(DG8). They share the same father with ‘Otisco’, ‘SX64’, 
but their overall yields were much lower, suggesting 
superior combining ability of ‘SV?’ for yield traits. This 
suggests that interspecific triploid hybridization in Salix 
can confer substantial gains in yield, and future breed- 
ing efforts will continue to exploit this potential. 

By correlating environment mean yields and IPC 
scores with edaphic and climatic variables from the test 
sites, we were able to demonstrate that increased tem- 
peratures were positively correlated with both mean 
yields and IPC1 scores, while the correlation between 
precipitation and mean yield was marginally significant. 
We have observed that the triploid hybrids tend to 
retain a leaf canopy much later in the season compared 
to other pedigrees. This could be an indication of lower 
sensitivity to frost or an ability to exploit solar resources 
at the tail end of the growing season. Labrecque & 
Teodorescu (2003) examined yields of two willow spe- 
cies at two contrasting sites in southern Québec and 
suggested that growth in that region may be limited by 
precipitation, but not likely growing season tempera- 
tures, based on reported optimal growing conditions in 
Sweden. Wang & Macfarlane (2012) analyzed poplar 
(Populus spp.) and shrub willow yields of cultivars 
sourced from Ontario and New York collections at two 
locations in Michigan and estimated that the difference 
in growing degree-days between northern and southern 
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locations could explain 60% of the variation in yield. In 
the UK, Aylott et al. (2008) analyzed yields of three wil- 
low GEN of differing pedigrees in relation to edaphic 
and climatic variables from 49 locations. In general, 
temporal rainfall patterns were found to be the signifi- 
cant factors affecting yields, but the response varied by 
GEN, with one being particularly sensitive to soil pH. 
We were unable to find significant correlations with soil 
pH or organic matter; however, because of the broad 
geographic range covered in this analysis, perhaps the 
greatest amount of variability among the test sites was 
dominated by climatic factors. The subject warrants 
additional investigation and testing in locations with 
higher temperatures and moderate rainfall may reveal a 
greater geographic production range than previously 
proposed (Walsh et al., 2003). 

GxE interactions have been reported elsewhere in 
recent studies examining willow yields in multiple loca- 
tions. In eastern Canada, Mosseler et al. (2014a) studied 
biomass traits of multiple clones from natural accessions 
of two species, S. discolor and S. eriocephala, planted on 
three sites of contrasting fertility. They found no signifi- 
cant rank changes between the two species groups, but 
large differences in genotypic sensitivities among the 
GEN within species, suggesting a high degree of 
intraspecific genetic variability in adaptation that could 
be exploited for biomass yield improvement. In Den- 
mark, Larsen et al. (2014) studied yields of Swedish and 
UK commercial cultivars across five locations and 
reported significant GxE interactions, but there was one 
clear winner, ‘Tordis’ (S. viminalis x S. schwerinii), with 
a mean yield of 6.7 Mg ha‘! yr`'. 
involving numerous GEN of three willow species across 
a large number of sites in Sweden, Ronnberg-Wastljung 
& Thorsén (1988) used a variant of the Finlay-Wilkinson 
regression to analyze GxE interactions. As in our analy- 
sis, linear regressions did not always adequately explain 
clonal response to ENV. The authors identified stable 
GEN using regression coefficients and mean yields, but 
they did not address specific adaptability. In poplar, 
Zalesny et al. (2009) performed an extensive analysis of 
53-79 GEN across four sites and multiple ages in the 
upper Midwestern USA. They assessed genotypic rank- 
ings and segregated GEN based on generalists (consis- 
tently high rankings) and specialist (especially variable 
rankings). More recent attempts to analyze GxE interac- 
tions in poplar have incorporated some of the more con- 
temporary statistical techniques used in plant breeding 
and agronomy. In an analysis of first-year growth of 
nine GEN in five locations in Spain, Sixto et al. (2011) 
applied the genotype and genotype x environment 
(GGE) biplot (Yan et al., 2000), which has some similari- 
ties with AMMI biplots and is used to visualize interac- 
tions (Gauch, 2006a; Gauch et al., 2008). More recently, 


In earlier work 


Sixto et al. (2014) analyzed biomass production at the 
end of a 3-year rotation for nine GEN in four of the 
locations from their previous study. They applied mixed 
model variants of stability parameter models including 
Finlay-Wilkinson and Eberhart—Russell and identified 
patterns of broad and specific adaptation among GEN 
and pedigrees. These mixed model approaches often 
consider ENV as random and have the advantage of 
handling unbalanced datasets (Smith et al., 2005). 

The importance of developing cultivars with specific 
adaptation can be emphasized when willow cultivation 
is placed into broader biological and geographical con- 
texts. In the UK and Swedish breeding programs, 
improved S. viminalis and hybrids such as S. vimi- 
nalis x S. schwerinii are highly productive and have 
demonstrated superior yields in numerous trials com- 
pared to other pedigrees (Aylott et al., 2008; Lindegaard 
et al., 2011; Larsen et al., 2014). In the USA, potato 
leafhopper (Empoasca fabae Harris) has proven to be 
extremely damaging to imported European S. viminalis 
cultivars (Smart & Cameron, 2008). Hybrid crosses of 
S. viminalis with S. miyabeana have introduced varying 
degrees of resistance to potato leafhopper, in addition to 
improved yields. We have observed more severe dam- 
age to ‘Tully Champion’ compared with ‘Fabius’ in NY 
locations (Gouker & Smart, 2015), which may have lim- 
ited its yield potential on sites of more southern latitude 
and more eastern longitude. Potentially lower potato 
leafhopper pressure in northerly ENV may explain the 
superiority of “Tully Champion’ and ‘SV1’ at those sites. 
In Canada, S. eriocephala and S. discolor have been evalu- 
ated as North American native species appropriate for 
biomass and phytoremediation applications (Mosseler 
et al., 2014a,b). However, testing in the USA has shown 
that natural accessions as well as improved cultivars of 
S. eriocephala can be highly susceptible to leaf rust 
(Melampsora spp) (Cameron et al., 2008; Serapiglia et al., 
2013). The improved S. eriocephala cultivar ‘S25’ in our 
analysis had the lowest overall yield, and it is possible 
that rust infection in combination with deer browse con- 
tributed to low performance. In Europe, Melampsora rust 
is a major growth-limiting pathogen for multiple willow 
species and species hybrids (Ahman, 1998; McCracken & 
Dawson, 2003; Aylott et al., 2008). Salix purpurea and 
hybrid S. koriyanagi x S. purpurea cultivars are also sus- 
ceptible to rust, but based on visual observations in the 
field and the yields of ‘Oneida’ and ‘Millbrook’ in this 
analysis, the S. purpurea x S. miyabeana diversity group 
likely has increased rust resistance. Breeding for rust 
resistance is an important focus in the USA given the 
European experience. 

Inclusion of some ENV in this analysis caused overall 
yields to be lower compared to other studies, but collec- 
tively they provided strong contrasts for discriminating 
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among GEN. In turn, this analysis likely provides a 
more realistic assessment of yield potentials on mar- 
ginal lands, given some of the exceptionally low-yield- 
ing ENV. As trials are often performed at experiment 
stations, where growing conditions are optimized or at 
least less variable, it may seem obvious that many 
reported yields are often on the high end and they may 
not be representative (Simmonds, 1991). This becomes 
even more relevant in the realm of dedicated bioenergy 
crops, as the premise is that these are purposed for mar- 
ginal lands; more specifically, where the economic 
returns are marginal for production of traditional field 
crops (Stoof et al., 2015), and thus, there is reduced com- 
petition between food and energy production. Evalua- 
tions of shrub willow yield potentials in the USA have 
been largely geographically restricted to regions where 
institutional knowledge exists, and the true extent of pro- 
duction potential has likely not be adequately tested 
(Walsh et al., 2003). Given the positive correlations of 
yield and IPC1 scores with factors relating to increased 
temperature, perhaps growing willow at lower latitudes 
will produce greater yields. But genotypic variation in 
water-use efficiency and drought resistance will be 
important and may help to inform breeding for improved 
adaptation to warmer climates (Bonosi et al., 2013). 

This study has demonstrated the importance of iden- 
tifying genotypic adaptability for developing cultivar 
recommendations and guiding future breeding efforts. 
Genotypes with high yields but varying sensitivities to 
environmental conditions have been identified as 
important check cultivars for testing the next generation 
of promising genotypes, which is currently underway. 
More work is needed in exploring the underlying envi- 
ronmental causes for the observed yields, and this will 
be the focus of a future analysis, which will incorporate 
an expanded dataset of North American willow yields. 
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